#!/bin/sh
####################################################################
#
#   Copyright 2003-2010 Peter Brommer, Daniel Schopf
#             Institute for Theoretical and Applied Physics
#             University of Stuttgart, D-70550 Stuttgart, Germany
#             http://www.itap.physik.uni-stuttgart.de/
#
#################################################################
#
#   This file is part of potfit.
#
#   potfit is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   potfit is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with potfit; if not, see <http://www.gnu.org/licenses/>.
#
#################################################################

wdir=`pwd`
list_types="0";
new_list="0";
f_arg="0";
s_arg="0";
e_file="";
saeng="0 0 0";
mycat="cat";
weight="1";
declare -a type_list;

# get command line options
while getopts 'c:e:flrs:w:?h' OPTION
do
    case $OPTION in
	c) c_string="$OPTARG";
            ;;
	e) e_file="$OPTARG";
            ;;
	l) list_types="1";
      	    ;;
	f) f_arg="1";
	    ;;
	s) s_arg="$OPTARG";
	    ;;
	r) scan_recursive="1";
	    ;;
	w) weight="$OPTARG";
	    ;;
	?) printf "\nUsage: %s: [-c list] [-e file] [-f] [-l] [-r] [-s list] <OUTCAR files>\n" $(basename $0) >&2
            printf "\n <OUTCAR files> is an optional list of files, if not given" >&2
	    printf "\n all files starting with OUTCAR will be scanned" >&2
	    printf "\n (it is possible to read gzipped files ending with .gz)\n" >&2
	    printf "\n -c <list>\t\tlist of chemical species to use\n" >&2
	    printf "\t\t\t(eg. \"vasp2force.sh -c Al=0,Mn=1,Pd=2\")\n" >&2
            printf " -e <file>\t\tspecify file to read single atom energies from\n" >&2
	    printf "\t\t\tif not found, \"0\" will be used for every atom type\n" >&2
	    printf " -f\t\t\tonly use the final configuration from OUTCAR\n" >&2
	    printf " -l\t\t\tlist all chemical species found in OUTCAR and exit\n" >&2
	    printf " -r\t\t\tscan recursively for OUTCAR files\n" >&2
	    printf " -s <list>\t\tcomma separated list of configurations to use\n" >&2
	    printf " -w <weight>\t\tchange weight for all configurations to <weight>\n" >&2
	    exit 2
	    ;;
    esac
done

# shift arguments by options and get the rest as $outcars
shift $(($OPTIND-1))
outcars="$*";

# if no outcars are given, search for them (recursively)
if [ "$outcars" == "" ]; then
    echo "No input files given on the command line. Searching for OUTCAR* ..." >&2;
    if [ "$scan_recursive" == "1" ]; then
    for i in `find . -name OUTCAR\* -print`; do
	outcars="${outcars} $i";
    done
    else
    for i in `find . -maxdepth 1 -name OUTCAR\* -print`; do
	outcars="${outcars} `basename $i`";
    done
    fi
    if [ "$outcars" != "" ]; then
	echo "Found the following files:$outcars" >&2;
    else
	echo "Could not find any OUTCAR files in this directory, aborting."
	break;
    fi
fi

for file in $outcars; do
    if [ ! -f $file ]; then
	continue;
    fi
    mycat="cat";
    zipped=`echo $file | awk 'sub(/\.gz$/,""){print $0}'`
    if [ "$zipped" != "" ]; then
	mycat="zcat";
    fi
    count=`$mycat $file | grep -c "TOTAL-FORCE"`;
    pr_conf="0";
    if [ "$f_arg" == "1" ]; then
	pr_conf="${pr_conf},$count";
    fi
    if [ "$s_arg" != "0" ]; then
	pr_conf="${pr_conf},$s_arg";
    fi
    types=`$mycat $file | grep VRHFIN | wc -l`;
    if [ $types != "0" ]; then
	name=(`$mycat $file | grep VRHFIN | awk '{ sub("=",""); sub(":",""); print $2; }'`);
    	if [ "$list_types" == "1" ]; then
	    echo "Found $types atom types in $count configurations in $file:";
	    for (( i=0; $i<$types; i++ )); do
		echo ${name[$i]} "= "$i;
	    done
	fi
    else
	types=`$mycat $file | grep TITEL | wc -l`;
	if [ $types != "0" ]; then
	    name=(`$mycat $file | grep TITEL | awk '{print $4; }'`);
    	    if [ "$list_types" == "1" ]; then
		echo "Found $types atom types in $count configurations in $file:";
		for (( i=0; $i<$types; i++ )); do
		    echo ${name[$i]} "= "$i;
		done
	    fi
	else
	    types=`$mycat $file | grep POTCAR | wc -l`;
	    types=$(($types/2));
		name=(`$mycat $file | grep POTCAR | awk '{print $3; }'`);
		if [ "$list_types" == "1" ]; then
		echo "Found $types atom types in $count configurations in $file:";
		for (( i=0; $i<$types; i++ )); do
		    echo ${name[$i]} "= "$i;
		done
	    fi
	fi
    fi
    if [ "$list_types" != "1" ]; then
	echo "There are $count configurations in $file" >&2
	if [ "$e_file" != "" ]; then
		if [ -f $e_file ]; then
		    saeng=`cat $e_file`;
		else
		    printf "Warning: $e_file could not be found - using 0 instead\n" >&2;
		fi
	fi

	if [ "X$pr_conf" == "X0" ]; then
	    pr_conf=1;
	    for (( i=2; $i<=$count; i++ )); do
		pr_conf="${pr_conf},$i";
	    done
	fi
	if [ ${#name[*]} != ${types} ]; then
		temp_c=${#name[*]};
		for (( i=${types};i<$temp_c;i++ )); do
			unset name[$i];
		done
	fi
	name_array=$( printf "%s," "${name[@]}" )
	$mycat $file | awk -v pr_conf="${pr_conf}" -v wdir="${wdir}" -v poscar="${poscar}" -v saeng="$saeng" -v name="$name_array" -v c_string="$c_string" -v file="$file" -v recursive="$scan_recursive" -v weight="$weight" '  BEGIN {
    OFMT="%11.7g"
#Select confs to print
    count=0;
    split(pr_conf,pr_arr,",");
    sub(/,$/,"",name);
    n_names=split(name,names,",")
    n_elements=split(c_string,el_names,",");
    sub(/=./,"",c_string);
    if (n_names>n_elements && c_string != "") {
	print "\nERROR:" > "/dev/stderr";
        print "There are more elements in the configuration than specified on the command line." > "/dev/stderr";
	print "Found "n_names" ("name") in "file" but expected only "n_elements" ("c_string")." > "/dev/stderr";
	exit;
    }
    if (c_string != "") {
	    asort(el_names);
	    for (i in el_names) {
		    split(el_names[i],temp,"=");
		    el_array[temp[1]]=temp[2];
		    el_inv_array[temp[2]]=temp[1];
	    }
	    for (i in el_names) {
	    el_number[names[i]]=el_array[names[i]];
    		}
    }
    for (i in pr_arr) pr_flag[pr_arr[i]]++;
    single_energy=0.;
    split(saeng,sae);
  };
  /ions per type/ {
    ntypes = split($0,a) - 4;
    for (i=1;i<=ntypes;i++) a[i]=a[i+4];
    for (i=1; i<=ntypes; i++) single_energy += a[i]*sae[i];
    for (i=2; i<=ntypes; i++) a[i]+=a[i-1];
    j=1; for (i=1; i<=a[ntypes]; i++) { if (i>a[j]) j++; b[i]=j-1; }
    }
# Marker for energy is two spaces between energy and without.
# Correct energy is energy(sigma->0)
  /energy  without/ {
    energy=($7-single_energy)/a[ntypes];
    delete stress;
  }
# Find box vectors
  /VOLUME and BASIS/ {
    for (i=1;i<=4;i++) getline;
    getline boxx; getline boxy; getline boxz;
    gsub("-"," -",boxx);
    gsub("-"," -",boxy);
    gsub("-"," -",boxz);
    split(boxx,boxx_v);
    split(boxy,boxy_v);
    split(boxz,boxz_v);
    scale=1.0;
}
  ($2=="kB") {
     for (i=1;i<=6;i++) stress[i]=$(i+2)/1602.;
}
  ($2=="TOTAL-FORCE") {
     count++;
     if (count in pr_flag) {
       print "#N",a[ntypes],1; #flag indicates whether to use forces or not
       printf "#C";
       if (c_string == "") {
       for (j=1;j<=ntypes;j++)
	       printf " %s",names[j];
       } else {
       for (j=1;j<=n_elements;j++) {
		printf " %s",el_inv_array[j-1];
	}
       }
       printf "\n";
       ("readlink -f "file) | getline x;

       print "## force file generated from file "x;

       printf "#X %13.8f %13.8f %13.8f\n",boxx_v[1]*scale,\
                   boxx_v[2]*scale,boxx_v[3]*scale;
       printf "#Y %13.8f %13.8f %13.8f\n",boxy_v[1]*scale,\
                   boxy_v[2]*scale,boxy_v[3]*scale;
       printf "#Z %13.8f %13.8f %13.8f\n",boxz_v[1]*scale,\
                   boxz_v[2]*scale,boxz_v[3]*scale;
       printf("#W %f\n",weight) ;
       printf("#E %.10f\n",energy) ;
       if ( 1 in stress )
         print "#S",stress[1],stress[2],stress[3],stress[4],stress[5],stress[6];
       print "#F";
     }
     getline; getline;
     for (i=1; i<=a[ntypes]; i++) {
	     if (c_string == "") {
       if (count in pr_flag)
	   printf("%d %11.7g %11.7g %11.7g %11.7g %11.7g %11.7g\n",
	       b[i],$1,$2,$3,$4,$5,$6);
       } else {
       if (count in pr_flag)
	   printf("%d %11.7g %11.7g %11.7g %11.7g %11.7g %11.7g\n",
	       el_number[names[b[i]+1]],$1,$2,$3,$4,$5,$6);
       }
     getline; }
  };'
    fi
done
